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ABSTRACT 

We present two-dimensional inviscid hydrodynamic simulations of a proto¬ 
planetary disk with an embedded planet, emphasizing the evolution of potential 
vorticity (the ratio of vorticity to density) and its dependence on numerical res¬ 
olutions. By analyzing the structure of spiral shocks made by the planet, we 
show that progressive changes of the potential vorticity caused by spiral shocks 
ultimately lead to the excitation of a secondary instability. We also demonstrate 
that very high numerical resolution is required to both follow the potential vor¬ 
ticity changes and identify the location where the secondary instability is hrst 
excited. Low-resolution results are shown to give the wrong location. We estab¬ 
lish the robustness of a secondary instability and its impact on the torque onto 
the planet. After the saturation of the instability, the disk shows large-scale non- 
axisymmetry, causing the torque on the planet to oscillate with large amplitude. 

The impact of the oscillating torque on the protoplanet’s migration remains to 
be investigated. 
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1. INTRODUCTION 

Well before extrasolar planets were discovered, Goldreich & Tremaine (1979; 1980) and 
Lin & Papaloizou (1986a,b) have speculated that tidal interactions between disks and embed¬ 
ded protoplanets would lead to planet migration. Ward (1997) suggested that two different 
types of migration could occur. One, called type I migration, is when the protoplanet mass 
is still small enough that the migration rate can be evaluated using linear analysis and is 
shown to be quite fast in the sense that migration timescale is shorter than the estimated 
buildup timescale of a giant protoplanet core. This has presented a serious problem for the 
formation of giant planets. The other, called type II migration, is when the protoplanet is 
massive enough that it opens gaps in the disk. The protoplanet is then locked into the disk’s 
viscous evolution timescale, which is typically long. 

The discovery of extrasolar planets with a large number of short-period planets very 
close to their parent stars seems to support a migration scenario since it is more likely 
that giant planet formation takes place at larger distances outside the snow line. This has 
generated renewed interests in the study of tidally induced planet migration (e.g., Kley 1999; 
Lubow, Seibert, & Artymowicz 1999; Nelson et ah 2000; D’Angelo, Kley, & Henning 2003). 
Both types of migration have been basically conhrmed by non-linear numerical simulations 
though most simulations to-date are performed with a prescribed viscosity. 

A recent study by Balmforth & Korycansky (2001) focused on the evolution of potential 
vorticity in the co-orbital region of a planet and they noted the possibility of secondary 
instabilities in the corotation resonance region, which lead to the formation of vortices. 
They showed that nonlinear dynamics of the co-orbital region have an important influence 
on the saturation of the corotation torque. 

In an earlier paper (Koller et ah 2003), we also investigated the nonlinear dynamics in 
the co-orbital region (within ~ 10 Roche lobe radii of the planet) through global nonlinear 
2D inviscid simulations. We found that the potential vorticity of the disk flow undergoes 
systematic changes, presumably caused by the spiral shocks generated by the planet though 
we did not give any detailed analysis of the shocks. We also found that the flow eventually 
becomes unstable due to the inflexion points in the potential vorticity prohle. Vortices are 
formed as a consequence. 

This paper is a close follow up to Koller et ah Our purpose is two fold. First, we present 
quantitative analyses of the ways that shocks change the potential vorticity of the flow. 
Second, through extensive numerical resolution studies, we conhrm the physical mechanism 
for exciting a secondary instability through inflexion points in the potential vorticity prohle. 
However, we now hnd that the instability is excited at a different location from what Koller 
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et al. presented. The phenomenon observed by Koller et ah is apparently due to limited 
numerical resolutions, hence artihcial. Nevertheless we are still able to demonstrate the 
existence of such a secondary instability, and its impact on the planet’s torque. We have 
also performed extensive numerical tests via different resolutions to check the validity of 
these new hndings. 

This paper is organized as follows. We first describe our numerical approach and how 
we set up our simulations in §2 and §3. We then present the analyses on shocks and how they 
change the potential vorticity in §4. The excitation of a secondary instability is described in 
§5. The influence of numerical resolution is presented in §6. Discussions are given in §7. 


2. NUMERICAL MODEL 

We assume that the protoplanetary disk is thin and can be described by the inviscid two- 
dimensional isothermal Euler equations in a cylindrical {r, 0} plane with vertically integrated 
quantities. The differential equations are the same as given in Kley (1998). Simulations are 
carried out mostly using a split dimension hydro code (Li & Li 2005) whose basic algorithm 
is based on the MUSCL-Hancock scheme (MRS) (Toro 1999). The methodology of MRS is 
extended here for the presence of source terms with two following main modihcations: the 
fluxes are computed using primitive instead of conserved variables, where van Leer flux lim¬ 
iter is applied to the slopes of the primitive variables, and an iterative approximate Riemann 
solver is used. Standard dimensional splitting as well as some source splitting are used. 
Radial sweeps with the angular flux derivatives omitted are alternated with angular sweeps 
with the radial flux derivatives omitted. All source terms except the planet’s gravitational 
force are included in the radial sweep. The planet force is applied in a separate sweep as 
the numerical gradient of the planet potential, so that the numerical curl of this gradient is 
zero, minimizing the potential vorticity contamination by the planet. 

We use the local co-moving angular sweep as proposed in the FARGO scheme of Masset 
(2000) and modihed in Li et al. (2001). The idea is basically to use a semi-Lagrangian form 
for the transport terms in the angular sweep. A constant velocity is subtracted from the 
angular transport velocity. This velocity is chosen to be close to the mean angular velocity, 
in order to move the data an integral number of angular cells in one time step. This is then 
corrected by an integral shift of the data. This enables us to use a much larger time step 
than would otherwise be possible. 

The MRS scheme requires two ghost cells in the radial sweep (the angular sweep is 
periodic). Rolding these ghost cells at the initial steady state produced the smallest boundary 



reflection of all boundary conditions tried. 


3. INITIAL SETUP 

The two-dimensional disk is modeled between 0.4 < r < 2. The planet is assumed to 
be on a hxed circular orbit at r = 1. A co-rotating frame is used and the positions of the 
central star and the planet are hxed at (r, 0) = (0,0) and (r^, ^p) = (l,7r) [acceleration due 
to frame rotation is also included, see Kley (1998)]. The mass ratio between the planet and 
the central star is /i = Mp/M*. Its Hill (Roche) radius is th = The disk is 

assumed as isothermal with a constant temperature throughout the simulation region. The 
isothermal sound speed, scaled by the Keplerian rotation speed at r = 1, is = H/r 
where H is the disk scale height. Time is measured by the orbital period at r = 1. 

We choose an initial surface density prohle with S(r) oc so that the ratio of 

vorticity to surface density, potential vorticity (PV = C); has a hat radial prohle C = (V x 
~ const. ~ 0.5. (Note that C, has a small deviation from being precisely 0.5 due to 
the hnite pressure gradient which slightly modihes from its Keplerian value.) With this 
initial condition, we avoid the generation of inhection points due to the rearrangement of C, 
distribution as it is carried by the stream lines (Balmforth & Korycansky 2001). 

We typically run the disk without the planet for 10 orbits hrst so that the disk is settled 

into a numerical equilibrium (very close to the analytic one). Then the planet’s gravitational 
potential is gradually “turned-on” over a 10-orbit period, allowing the disk to respond to the 
planet potential gradually. Furthermore, the planet’s potential is softened by an approximate 
three-dimensional treatment. The basic idea is that, at any location {r, 0}, the matter with 
surface density S is distributed (divided into many cells) vertically over many scale heights 
(say ±6iif) according to a Gaussian isothermal prohle. Radial and azimuthal forces exerted 
by the planet at a specihc {r, 0, z} are calculated in this 3D fashion. To get the 2D radial 
and azimuthal forces for a specihc {r, 0}, one then integrates all the cells along ; 2 . Comparing 
this pseudo-3D treatment with the usual 2D treatment where a hxed smoothing distance is 
used, we hud that it reduces the force close to the planet (within ~ 2r/f) by about a factor of 
2 but converges to the usual 2D force beyond ~ 5rH- Although there is some accumulation 
of gas near the planet, we do not allow any disk gas to be accreted onto it. Self-gravity of the 
disk is not included. Runs are made using several radial and azimuthal grids to study the 
inhuence of resolution, which ranges from {ur x n^) = 200 x 800 to 1200 x 4800. Simulations 
typically last several hundred orbits at r = 1. 

One key feature of our simulations is that they are performed at the inviscid limit, i.e.. 
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we do not explicitly include a viscosity term, though numerical viscosity is inevitable and is 
needed to handle shocks. In addition, our simulations tend to be of higher resolution. With 
a radial Ur = 800 (some runs go up to = 1200, = 4800) and 0.4 < r < 2, the diameter 

of the Hill sphere of a planet with p = 10“^ is resolved by ~ 32 cells in each direction. 

4. SHOCK STRUCTURE and PV CHANGES 

We will use one run in this section to facilitate discussion. This run has p = 10“^, 
Cs = 0.05, th = 0.0322, and a x = 800 x 3200 resolution. Figure 1 shows the surface 
density (multiplied by so that the background surface density is unity initially) in the 
{r, 0} plane at f = 100 orbits. It shows the characteristic density enhancement at the planet 
and at the two spiral shocks. Two density “depressions” at |Ar| ~ 3 — 4:rH are developed as 
angular momentum is gradually deposited there. Two density bumps (enhancements) occur 
broadly over |Ar| ~ 5 — 8rH- Here, Ar = r — Vp. Some streamlines (in the co-rotating frame) 
are plotted as well. 


4.1. PV Change by Shocks 

Figure 2 shows the azimuthally averaged {() prohle from t = 40,100,160, 220, and 300 
orbits. (The initial value ~ 0.5 is subtracted.) Several features can be seen. First, there are 
very small changes in {() for |Ar| < 2rH- Second, at larger distances from the planet, {() 
prohle changes progressively with time and shows both increase and decrease from its initial 
value. An interesting feature is that even though both “peaks” and “valleys” are growing 
with time, two locations around Ar = ±4.5rH have no change in {(). A natural question is 
what causes these features. 

The how around the planet can be separated into several sub-regions depending on 
their streamline behavior (e.g., see Masset 2002), including the horseshoe (or librating) 
region with |Ar| < th, the separatrix region th < |Ar| < vT2r//, and the streaming region 
|Ar| > y/VlrH- Typically, spiral shocks do not occur until the relative velocity between the 
disk how and the planet becomes larger than the sound speed. This roughly translates to 
|Ar| ^ H. So, a critical parameter is H/th ~ Cg/rn, which gives an indication whether 
the spiral shock is starting from within the planet’s Roche lobe or not (see also Korycanski 
& Papaloizou 1996). In the limit of small Cg/rn, how around the planet becomes quite 
complicated (e.g., see Tanigawa & Watanabe 2002). Another added yet potentially very 
important complication is that any numerical error coming out nearing the planet might be 
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strongly amplified by these shocks, giving numerical artifacts which can be physically unreal. 
So, we have mostly used a relatively higher sound speed of Cg = 0.05 (> r^) in this study. 
Higher sound speeds mean weaker shocks, which typically also mean slower changes in the 
disk properties. 

As discussed in Koller et al. (2003), spiral shocks are capable of destroying the conser¬ 
vation of PV. For the parameters used in this run, we do not expect shocks within ±2r^,, 
which is conhrmed by the fact that {() stays roughly at the initial value. Once in the shocked 
region, we can use a very useful formula for calculating the PV jump across a (curved) shock, 
which is given in Kevlahan (1997), after some manipulation, 

^ Ml dr E 2 ’ ^ 

where Mj_ is the perpendicular Mach number {u±/cs) which is equal to (S 2 /Si)^/^ for isother¬ 
mal shocks. Si and S 2 are the pre- and post-shock surface densities. The direction r is 
tangential to the shock front and is defined to be pointing away from the planet. According 
to Eq. (1), the sign of A( is determined only by the gradient of M±_ along the shock front. 

Figure 3 shows M±_ as a function of radial distance to the planet (only one side is 
shown). This is obtained by Erst identifying the high compression region through calculating 
the velocity compression. After finding the flow impact angle to that surface, one can then 
calculate M±_ at the shock. Differentiating along that surface gives dM^/dr. This figure 
shows the general behavior of what one might expect of a spiral shock produced by the 
planet. As one moves radially away from the planet, the flow goes from a non-shocked state 
to having shocks, which means that M_l gradually increases from being < 1 to 1 when the 
shock actually starts. Very far away from the planet, the pitch angle between the spiral 
wave and the background flow becomes small enough that M_i_ is expected to become less 
than one again. So, as a function of radial distance from the planet, M\_ is expected to 
increase first, reaching its peak, then starting to decrease, it eventually drops to be smaller 
than one at some large distance. Consequently, using eq. (1), we can expect that should 
be positive first, becomes zero when M_i_ reaching its peak, then is negative when M_i_ is 
decreasing. This matches quite well with the behavior of {Q given in Fig. 2. Note that over 
the whole evolution, shocks are quite steady with very slow evolution, indicated by the small 
“outward” movement of the inflexion points (where AC, ~ 0). All these features are quite 
generic and we have confirmed this by using several different planet masses. 

From Fig. 3, we can see that M_i_ peaks at Ar —S.Sth, but the point where AC = 0 
is at —AAth based on Fig. 2. To quantitatively match the two profiles in terms of radial 
distances to the planet, one has to take into the account the fact that the radial location 
where shock occurs (and the flow’s PV is changed) is different from where these PV changes 
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are eventually deposited. This is illustrated in Fig. 1. Correcting for this systematic shift 
in radial locations, the agreement between the azimuthally averaged {() prohle and the M± 
prohle is quite good. 

To summarize, the spiral shocks emanating close to the planet cause systematic changes 
to the PV distribution of the disk flow. As the disk rotates, the disk flow repeatedly passes 
through these shocks so that changes in PV then accumulate. The sign of PV change is 
determined by the shock structure, and the amplitude of change is determined both by the 
shock structure [eq. (1)] and how many times the flow passes through the shocks in a given 
time [~ — l)At/27r]. We have identihed the physical cause for having both positive 

and negative changes in PV, in addition to the existence of an inflexion point. We defer the 
quantitative comparison of the amplitude of PV changes between simulations and eq. (1) to 
a future study. 


5. SECONDARY INSTABILITY 
5.1. Instability and Vortices 

Having established the prohle of we now turn to the issue of secondary insta¬ 

bility. This type of instability is usually associated with the existence of inhexion points 
(extrema) in the PV prohle (see Drazin & Reid 1981) and often requires satisfying a thresh¬ 
old condition. As shown in Fig. 2, the PV prohle has extrema points, and the magnitude of 
the slope between the peak at Ar 3.5rH and the valley at Ar ^ 4.6r// is increasing with 
time (similarly on the other side of the planet). It is then expected that these changes will 
become large enough that a threshold is reached and secondary instabilities will be excited. 

To demonstrate this, we have made runs out to several hundred orbits (400 — 800 orbits 
depending on resolution). Figure 4 shows the evolution of PV in the {r, 0} plane for a run 
with fi = 10“"^, Cs = 0.05, and Ur x = 400 x 1600. Indeed, a secondary instability is 
observed to occur in the outer “valley” of the PV prohle at t ~ 180 orbits. (The other 
valley in the ~ region also becomes unstable at a later time.) The negative PV 

regions correspond to the two density bumps, which are roughly axisymmetric early on. 
When the instability is excited, density bumps break up into ~ 5 anticyclonic vortices (even 
higher density blobs). Over the next tens of orbits, as the instability evolves and saturates, 
these density blobs merge, forming density enhancements that are non-axisymmetric on large 
scales. This is illustrated in Fig. 5 at t = 500 orbits. Note that density depression regions 
also become non-axisymmetric. 



5.2. Torque Evolution 


As discussed in Koller et al. (2003), one primary goal of studying secondary instability 
is to understand its impact on the torque evolution of disk gas on the planet. Figure 6 shows 
the torque history on the planet for a run with the same parameters as in Figure 4. The rapid 
oscillations in the torque history starting at t ~ 180 indicates the initiation of the instability 
(consistent with Fig. 4). As the higher density blobs pass by the planet, they exert stronger 
torques. Furthermore, Figure 6 shows that the oscillation amplitude grows with time, which 
is caused by the gradual density increase in the PV “valleys” over long timescales. To better 
understand the origin of fast oscillations, we have calculated the torques separately from 
three different parts of the disk, which are inner: Ar < Sr^, middle: |Ar| < Sr^, and outer: 
Ar > Sth, respectively. Figure 7 shows their individual contributions over a time interval 
from t = 490 — 512. We can see that both inner and outer regions have a periodic large 
amplitude variation. These are caused by the orbital motion of the non-axisymmetric disk 
flow at Ar ~ —5rH and Ar ~ 7.5rH, respectively (see Fig. 5). In fact, the torque oscillation 
periods from these two regions are precisely the expected orbital periods at those radii. Since 
these two periods are different, their sum gives an erratic appearance (Fig. 6). 

One might speculate that even with these large amplitudes and rapid oscillations, there 
seems to be a mean value that is negative and has a small amplitude, consistent with the 
type I migration expectations. But since our simulations have the planet moving on a hxed 
circular orbit, it is premature to conclude that these oscillations will not change the type 
I migration picture. Studies allowing the planet to migrate under the influence of these 
oscillating torques is under way to address this important issue. 

Here, we have mostly concentrated on demonstrating the existence of such a secondary 
instability and its potential impact on the planet’s torque. We have not quantitatively 
analyzed what is the threshold condition in the PV prohle that excites this instability. Such 
a threshold obviously depends on the planet mass, shock structure, sound speed of disk gas, 
and viscosity in the disk, etc. We have assumed the inviscid limit. Sufficient viscosity could 
potentially remove the buildup of the PV’s peaks and valleys made by shocks, hence never 
allowing the prohle to reach the critical destabilization level. 


6. RESOLUTION STUDY 

The results we presented here, especially Figs. 2 and 4, are different from the hnd- 
ings in Koller et al. (2003, see their Fig. 3, also curve A in our Fig. 8). Results there 
showed additional “dips” interior (i.e., closer to the planet) to the main positive peaks in 



- 9 - 

{(). Furthermore, those “dips” were shown to deepen with time and eventually became 
unstable. 

Figure 8 shows a comparison of runs with 4 different resolutions but the same set of 
physical parameters. Even though different curves of {() converge at large Ar (> Gr^) 
because the shock is weak and the gradient is small there, the behavior at small Ar (where 
the shock is strong) shows large differences. Comparing curves A and C, for example, we 
can see that when the resolution is low and shocks are not well resolved, a larger part of the 
disk was affected (see the region of Ar ~ 0.5 — 2.5rH in curve A). We were able to reproduce 
the results of Koller et al. in onr low-resolution runs bnt our high-resolution studies indicate 
that the “dips” observed in Koller et al. are getting narrower and shallower with higher 
resolutions. This means that the previous results showing the inner “dips” are nnmerical 
artifacts, not true physical effects. Note that the PV within the planet’s Roche lobe is not 
conserved dne to both the switch-on of the planet’s mass (no matter how slowly) and the 
nnmerical error in implementing the planet’s potential. This is indicated by Fig. 9, where 
PV in the {r, 0} plane at t = 100 is shown for two different resolutions. Small PV changes 
are visible coming ont of the planet region. For the low resolntion case, the PV change from 
the planet is “spread” over a wide region. It is then amplified by poorly resolved shocks. 
Such amplihcation eventually leads to an instability, nnfortnnately, in mnch the same spirit 
we have discussed here. For the high resolution run, the shock and the PV change from the 
planet is well separated. The erroneous amplihcation does not occur. 

It is seen that very high resolutions (800 x 3200 to 1200 x 4800) are needed to obtain 
convergence in the {() prohle (curves C and D). We emphasize that even though there are still 
some minor differences for the high resolution runs we presented here, the overall feature 
of having a positive peak and a negative valley is qnite consistent among all resolutions. 
Fnrthermore, the instability now comes out of the valley region, which typically has a mnch 
better convergence than other parts of the ( prohle. In addition, we have also made a run 
using pure Lax-Wendroh scheme with rirXn^ = 800 x 3200 (curve E). Comparing with curve 
C, it gives very similar results. It also shows the instability at later times (not shown here), 
though the MHS scheme gives sharper and smoother shocks than hybrid schemes. 


7. CONCLUSIONS 

We have carried out high-resolution two-dimensional hydrodynamic disk simulations 
with one embedded protoplanet. We hnd that the total torqne on the planet, caused by 
tidal interactions between the disk and the planet, can be divided into two stages. First, it 
is negative, and slowly varying, consistent with the type I migration expectation. Second, 
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it shows large amplitude and very fast oscillations, due to the excitation of an instability 
which hrst breaks up the axisymmetric density enhancement into higher density blobs, e.g., 
vortices. These vortices then merge, forming large-scale non-axisymmetric density enhance¬ 
ments. This non-axisymmetry causes the torque to continuously oscillate. 

In Koller et ah (2003), we were misled by a spurious feature in the PV prohle, due to an 
inadequate numerical resolution, which eventually became unstable. Despite the inaccurate 
location where a secondary instability is initiated, the physical explanation proposed by 
Koller et al. for exciting such an instability is well founded. We now have a self-consistent 
picture from performing very high resolution simulations to understand the shock structure 
and the PV prohle. Further studies allowing the planet to migrate under the inhuence of 
these torques will be necessary and very interesting. 


This research was performed under the auspices of the Department of Energy. It was 
supported by the Laboratory Directed Research and Development Program at Los Alamos 
and by LANL/IGPP. 
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Fig. 1.— Surface density, multiplied by in the {r, 0} plane at t = 100 orbits for a run 
with /i = 10“^, Cs = 0.05, and Ur xn^ = 800 x 3200. (The initial background surface density 
is unity, after multiplying r^'®.) The planet is located at Ar = 0 and 0 = tt. Two spiral 
shocks cut through the disk, bending stream lines. The PV changes induced by the shock 
are transported to different radii from their prodnction site. 
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Fig. 2.— Azimuthally averaged {() profile at t = 40 (solid curve), 100,160,220, and 
300 (dot — dashed curve) orbits for a run with the same parameters as in Fig. 1. The 
planet is located at Ar = 0. The initial radial profile ( ~ 0.5 has been subtracted. Note the 
progressive growth of “peaks” and “valleys” in the {() prohle. Also, {() in the planet region 
remains roughly unchanged. 



Fig. 3.— Perpendicular Mach number M± along the spiral shock front for the run shown in 
Fig. 1 at t = 100 orbits. Note that the radial distances here mark the PV production sites, 
to which systematic shifts shown in Fig. 1 have to be added when compared with the radial 
locations in Fig. 2. 
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Fig. 4.— PV in the {r, 0} plane at t = 160,180,200, and 240 orbits (panels A-D) with 
rir X = 400 x 1600. (Results are similar for other resolutions.) Panels B and C show 
the initiation of a secondary instability in the “valley” regions with |Ar| = 5 — 7rH- This 
instability breaks the flow into higher density blobs (“vortices”) which are indicated by 
the isolated PV depressions (panels B and C). These vortices merge at late times, forming 
large-scale asymmetries in the azimuthal direction (panel D) . 


Fig. 5.— Surface density, multiplied by in the {r, 0} plane at t = 500 orbits. (Initial 
surface density is unity at r = 1, after multiplying r^'®.) After the instability has saturated, 
the density distribution becomes strongly non-axisymmetric. When they move around the 
disk, the torque on the planet oscillates. 



Orbits 


Fig. 6.— The total torque on the planet vs. time from a run having the same parameters 
in Fig. 4. Note that once the flow becomes unstable after t = 180, the total torque shows 
large amplitude oscillations. 
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Fig. 7.— A zoomed-in and “decomposed” view of Fig. 6 during orbits t = 490 — 512. The 
solid, dotted and dashed curves are torques from three different parts of the disk: Ar < 3r//, 
|Ar| < 3r//, and Ar > 3r//, respectively. The periodic oscillations in torque prohles (solid 
and dashed curves) are caused by a disk with non-axisymmetric density prohles rotating 
around the planet (see Fig. 5). When the higher density region moves close to the planet, 
it gives a stronger interaction. The oscillation periods in torques are the same as their 
respective orbital periods. 
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Fig. 8. — Azimuthally averaged (C) prohle at t = 160 orbits with yU = 10“^ and Cg = 0.05, 
but different resolutions: curves (A-D) use the MHS scheme with Ur xn^ = 200 x 800,400 x 
1600, 800 X 3200,1200 x 4800, respectively. Curve E uses the pure Lax-Wendroff scheme 
with Ur X = 800 x 3200. Curves C and D are quite similar to each other, indicating 
convergence. 
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Fig. 9.— PV in the {r, 0} plane at t = 100 orbits with rir x = 200 x 800 (left) and 
800 X 3200 (right). The planet is located at Ar = 0 and (j) = n. Only one side is shown. 
Resolving the shock is important in determining the starting point of the shock and where 
the PV changes will be deposited. 
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